* Create results reported in Section 6.2 - Scenario 1 (replace bottom 10% of doctors)

clear all

use "$savedata/masterdata.dta", replace

gen nstemi3 = ami3-stemi3

keep if sample25==1
gen vol = vol25

gen ps2 = 2
replace ps2 = 1 if stemi2==1 
table ps2, c(mean all_death30 mean all_death365)

foreach x in 1 2{
gen one`x' = 0
replace one`x' = 1 if ps2==`x'

egen n_efe_n`x' = sum(one`x'), by(doctor_id)
egen vol`x' = sum(one`x'), by(doctor_id)
}

foreach x in 1 2{
gen one`x'_2017 = 0 
replace one`x'_2017 = 1 if ps2==`x' & finyear==2017
}


foreach x in 1 2{
reghdfe survive365 c.prevyear_cost if ps2==`x', absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual`x', residuals

xtreg residual`x' c.stemi3 c.nstemi3 c.std_nonami3, fe i(doctor_id)
predict docfe_n`x', u
gen sigma_u_n`x' = e(sigma_u)
gen sigma_e_n`x' = e(sigma_e)
}

foreach x in 1 2{
gen signal_n`x' = ((sigma_u_n`x'^2) / ((sigma_u_n`x'^2) + ((sigma_e_n`x'^2)/n_efe_n`x')))
gen adj_docfe_n`x' = docfe_n`x'*signal_n`x'
}


* Also create the standard measure
reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step


* Just look at 2017
keep if finyear==2017

preserve
* 1) Replace bottom 10% of doctors and calculate mortality "gain"
collapse (mean) *docfe* (sum) vol_2017 = one, by(doctor_id)

xtile rank = adj_docfe, n(1069)
replace rank = 1070-rank

xtile deciles = adj_docfe, n(10)

gen docs_drop = 0
replace docs_drop=1 if deciles==1

gen pats_drop = vol_2017*docs_drop

gen extra_survival = pats_drop*(-adj_docfe)

egen gain = sum(extra_survival)
su gain /* This gives you the reduction in deaths reported in Section 6.2 */

restore

preserve
* 2) Reallocate doctors towards their comparative advantage (regardless of the hospital they work in)
collapse (mean) docfe* sigma* n_efe* (sum) one1 one2 one1_2017 one2_2017, by(doctor_id)

foreach x in 1 2{
gen signal_n`x' = ((sigma_u_n`x'^2) / ((sigma_u_n`x'^2) + ((sigma_e_n`x'^2)/n_efe_n`x')))
gen adj_docfe_n`x' = docfe_n`x'*signal_n`x'
su adj_docfe_n`x', det
}

* Now sort doctors by the difference (assume that people with missing values for low types have a zero difference)
gen docfe_high = adj_docfe_n1
gen docfe_low = adj_docfe_n2
egen temp = mean(adj_docfe_n2) 
replace docfe_low = temp if docfe_low==.
drop temp
egen temp = mean(adj_docfe_n1) 
replace docfe_high = temp if docfe_high==.
drop temp

gen docfe_diff = docfe_high - docfe_low

* Sort by their comparative advantage
sort docfe_diff

xtile rank = docfe_diff, n(1069)
replace rank = 1070 - rank

* Count capacity; total and by type of patient
gen total17 = one1_2017 + one2_2017
gen high_vol17 = one1_2017
gen low_vol17 = one2_2017

keep rank docfe_low docfe_high docfe_diff low_vol high_vol total doctor_id

sort rank

* Count total number of patients treated by each type, and cumulatively, up to each rank
gen add_count = total if rank==1

forval x=2(1)1184{

	local i=`x'-1
	
		qui gen temp = add_count if rank==`i'
		qui egen temp_count = max(temp)
		qui replace add_count = temp_count + total if rank==`x'
		qui drop temp 
		qui drop temp_count
		
	}	
	
egen total_high_type = sum(high_vol)

gen high_assign = 0
replace high_assign = 1 if add_count <= total_high_type


* Count the number of high and low cases they now have
gen new_high = total17*high_assign
gen new_low = total17 - new_high

gen change_high = new_high - high_vol
gen change_low = new_low - low_vol

* Calculate extra survival

gen extra_survival = change_high*docfe_diff
egen gain = sum(extra_survival)

su gain /* Gives you lives saved as a result of reallocation exercise 2 (unconstrained) */

save "$savedata/gains_unconstrained.dta", replace

restore

* 3) Constrained reallocation (within trust)

** First I need to create a record of the hospitals that doctors work in in 2017, and the volume of patient types in those hospitals
preserve
keep if finyear==2017
collapse (sum) one one1_2017 one2_2017, by(doctor_id trust_code)
rename one total_wt
rename one1_2017 total1_wt
rename one2_2017 total2_wt

save "$savedata/doc_hosps_2017.dta", replace

restore 

* Then use this to reallocate doctors

preserve
collapse (mean) adj* (sum) one1 one2 one1_2017 one2_2017, by(doctor_id)

gen diff_adj_docfe = adj_docfe_n1 - adj_docfe_n2

* Now sort doctors by the difference (assume that people with missing values for low types have a zero difference?)

gen docfe_high = adj_docfe_n1
gen docfe_low = adj_docfe_n2
egen temp = mean(adj_docfe_n2) 
replace docfe_low = temp if docfe_low==.
gen docfe_diff = docfe_high - docfe_low

* Sort by their comparative advantage
sort docfe_diff

* Count total capacity
gen total = one1 + one2
gen high_vol = one1
gen low_vol = one2

gen total17 = one1_2017 + one2_2017
gen high_vol17 = one1_2017
gen low_vol17 = one2_2017

keep if total17>0

xtile rank = docfe_diff, n(1069)
replace rank = 1070 - rank

keep rank docfe_low docfe_high docfe_diff low_vol17 high_vol17 total17 doctor_id


* Now match in the hospitals - we want an observation for each doc-hospital combination
merge 1:m doctor_id using "$savedata/doc_hosps_2017.dta"

sort trust_code rank

bys trust_code (rank): gen trust_count = sum(total_wt)

rename total_wt total17_wt
rename total1_wt high_vol17_wt
rename total2_wt low_vol17_wt

egen total_high_type_wt = sum(high_vol17_wt), by(trust_code)
egen total_low_type_wt = sum(low_vol17_wt), by(trust_code)

gen high_assign_wt = 0
replace high_assign_wt = 1 if trust_count <= total_high_type_wt

* Count the number of high and low cases they now have
gen new_high = total17_wt*high_assign
gen new_low = total17_wt - new_high

gen change_high = new_high - high_vol17_wt
gen change_low = new_low - low_vol17_wt

* Calculate extra survival

gen extra_survival = change_high*docfe_diff

egen gain = sum(extra_survival) /* Gives you lives saved from within-trust reallocation */

su gain

save "$savedata/gains_constrained.dta", replace

restore
